R语言生存分析的实现
💡专注R语言在🩺生物医学中的使用
欢迎加入🐧QQ交流群:613637742
生存分析是临床常用统计方法,一旦和时间扯上关系,分析就变得复杂多了,此时不再是单一的因变量,还需要考虑时间给因变量和自变量带来的各种影响。
本次主要演示R语言做生存分析的一些方法。比如寿命表、K-M曲线、logrank检验。后续还会给大家介绍Cox回归、时依系数和时依协变量的Cox回归、生存曲线的可视化等内容。
本推文不涉及理论,只有实操,想要了解生存分析的理论的请自行学习。不涉及理论,并不代表理论不重要,在以后的机器学习和临床预测模型的相关推文中,会经常用到这些理论,建议大家学习一下。
生存过程的描述
library(survival)
library(survminer)
使用survival
包中的lung
数据集用于演示,这是一份关于肺癌患者的生存数据。time
是生存时间,以天为单位,status
是生存状态,1代表删失,2代表死亡。但是一般在生存分析中我们喜欢用1代表死亡,用0代表删失,所以我们更改一下(其实不改也可以,你记住就行)。
df <- lung
df$status <- ifelse(df$status == 2,1,0)
str(df)
## 'data.frame': 228 obs. of 10 variables:
## $ inst : num 3 3 3 5 1 12 7 11 1 7 ...
## $ time : num 306 455 1010 210 883 ...
## $ status : num 1 1 0 1 1 0 1 1 1 1 ...
## $ age : num 74 68 56 57 60 74 68 71 53 61 ...
## $ sex : num 1 1 1 1 1 1 2 2 1 1 ...
## $ ph.ecog : num 1 0 0 1 0 1 2 2 1 2 ...
## $ ph.karno : num 90 90 90 90 100 50 70 60 70 70 ...
## $ pat.karno: num 100 90 90 60 90 80 60 80 80 70 ...
## $ meal.cal : num 1175 1225 NA 1150 NA ...
## $ wt.loss : num NA 15 15 11 0 0 10 1 16 34 ...
首先把生存时间和生存状态用Surv()
放到一起,可以看到有+
的就是截尾数据。
Surv(time = lung$time, event = lung$status)
## [1] 306 455 1010+ 210 883 1022+ 310 361 218 166 170 654
## [13] 728 71 567 144 613 707 61 88 301 81 624 371
## [25] 394 520 574 118 390 12 473 26 533 107 53 122
## [37] 814 965+ 93 731 460 153 433 145 583 95 303 519
## [49] 643 765 735 189 53 246 689 65 5 132 687 345
## [61] 444 223 175 60 163 65 208 821+ 428 230 840+ 305
## [73] 11 132 226 426 705 363 11 176 791 95 196+ 167
## [85] 806+ 284 641 147 740+ 163 655 239 88 245 588+ 30
## [97] 179 310 477 166 559+ 450 364 107 177 156 529+ 11
## [109] 429 351 15 181 283 201 524 13 212 524 288 363
## [121] 442 199 550 54 558 207 92 60 551+ 543+ 293 202
## [133] 353 511+ 267 511+ 371 387 457 337 201 404+ 222 62
## [145] 458+ 356+ 353 163 31 340 229 444+ 315+ 182 156 329
## [157] 364+ 291 179 376+ 384+ 268 292+ 142 413+ 266+ 194 320
## [169] 181 285 301+ 348 197 382+ 303+ 296+ 180 186 145 269+
## [181] 300+ 284+ 350 272+ 292+ 332+ 285 259+ 110 286 270 81
## [193] 131 225+ 269 225+ 243+ 279+ 276+ 135 79 59 240+ 202+
## [205] 235+ 105 224+ 239 237+ 173+ 252+ 221+ 185+ 92+ 13 222+
## [217] 192+ 183 211+ 175+ 197+ 203+ 116 188+ 191+ 105+ 174+ 177+
如果只是想要描述一下这份生存数据,可以使用寿命表法或者K-M曲线,在R中可以通过survfit()
实现。
# 构建生存曲线
fit <- survfit(Surv(time, status) ~ 1, data = df)
# 寿命表,surv_summary比默认的summary()更好
surv_summary(fit)
## time n.risk n.event n.censor surv std.err upper lower
## 1 5 228 1 0 0.99561404 0.004395615 1.0000000 0.98707342
## 2 11 227 3 0 0.98245614 0.008849904 0.9996460 0.96556190
## 3 12 224 1 0 0.97807018 0.009916654 0.9972662 0.95924368
## 4 13 223 2 0 0.96929825 0.011786516 0.9919508 0.94716300
## 5 15 221 1 0 0.96491228 0.012628921 0.9890941 0.94132171
## 6 26 220 1 0 0.96052632 0.013425540 0.9861367 0.93558107
## 7 30 219 1 0 0.95614035 0.014184183 0.9830945 0.92992527
## 8 31 218 1 0 0.95175439 0.014910735 0.9799794 0.92434234
## 9 53 217 2 0 0.94298246 0.016284897 0.9735659 0.91335978
## 10 54 215 1 0 0.93859649 0.016939076 0.9702809 0.90794671
## 11 59 214 1 0 0.93421053 0.017574720 0.9669508 0.90257880
## 12 60 213 2 0 0.92543860 0.018798180 0.9601711 0.89196244
## 13 61 211 1 0 0.92105263 0.019389168 0.9567281 0.88670745
## 14 62 210 1 0 0.91666667 0.019968077 0.9532533 0.88148430
## 15 65 209 2 0 0.90789474 0.021093908 0.9462168 0.87112471
## 16 71 207 1 0 0.90350877 0.021642644 0.9426590 0.86598451
## 17 79 206 1 0 0.89912281 0.022182963 0.9390770 0.86086855
## 18 81 205 2 0 0.89035088 0.023240987 0.9318456 0.85070391
## 19 88 203 2 0 0.88157895 0.024272607 0.9245323 0.84062118
## 20 92 201 1 1 0.87719298 0.024779731 0.9208475 0.83560803
## 省略一部分结果。。。
## 178 791 9 1 0 0.07831533 0.314346559 0.1450170 0.04229358
## 179 806 8 0 1 0.07831533 0.314346559 0.1450170 0.04229358
## 180 814 7 1 0 0.06712742 0.350176075 0.1333431 0.03379322
## 181 821 6 0 1 0.06712742 0.350176075 0.1333431 0.03379322
## 182 840 5 0 1 0.06712742 0.350176075 0.1333431 0.03379322
## 183 883 4 1 0 0.05034557 0.453824434 0.1225342 0.02068546
## 184 965 3 0 1 0.05034557 0.453824434 0.1225342 0.02068546
## 185 1010 2 0 1 0.05034557 0.453824434 0.1225342 0.02068546
## 186 1022 1 0 1 0.05034557 0.453824434 0.1225342 0.02068546
画出生存曲线,横坐标是生存时间,纵坐标是生存率。
ggsurvplot(fit,
conf.int = TRUE, # 可信区间
palette= 'blue', # 更改配色
surv.median.line = "hv", # 中位生存时间
ggtheme = theme_bw() # 更改主题
)
生存过程的比较
如果通过某个变量把数据分为多组,然后检验不同组别之间的生存时间(生存曲线)有无差别,则可以通过logrank检验或者breslow检验。
在R语言中通过survdiff()
实现logrank检验。
fit <- survdiff(Surv(time, status) ~ sex, data = df)
fit
## Call:
## survdiff(formula = Surv(time, status) ~ sex, data = df)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=1 138 112 91.6 4.55 10.3
## sex=2 90 53 73.4 5.68 10.3
##
## Chisq= 10.3 on 1 degrees of freedom, p= 0.001
可以用神包broom
提取数据:
broom::tidy(fit)
## # A tibble: 2 × 4
## sex N obs exp
## <chr> <dbl> <dbl> <dbl>
## 1 1 138 112 91.6
## 2 2 90 53 73.4
broom::glance(fit)
## # A tibble: 1 × 3
## statistic df p.value
## <dbl> <dbl> <dbl>
## 1 10.3 1 0.00131
对于不同组别之间生存曲线的检验,也可以通过K-M图示的方法:
fit.logrank <- survfit(Surv(time, status) ~ sex, data = df)
surv_summary(fit.logrank) # 可以查看寿命表
## time n.risk n.event n.censor surv std.err upper lower
## 1 11 138 3 0 0.97826087 0.01268978 1.0000000 0.95423012
## 2 12 135 1 0 0.97101449 0.01470747 0.9994124 0.94342350
## 3 13 134 2 0 0.95652174 0.01814885 0.9911586 0.92309525
## 4 15 132 1 0 0.94927536 0.01967768 0.9866017 0.91336116
## 5 26 131 1 0 0.94202899 0.02111708 0.9818365 0.90383547
## 6 30 130 1 0 0.93478261 0.02248469 0.9768989 0.89448204
## 7 31 129 1 0 0.92753623 0.02379334 0.9718155 0.88527450
## 8 53 128 2 0 0.91304348 0.02627035 0.9612864 0.86722163
## 9 54 126 1 0 0.90579710 0.02745220 0.9558688 0.85834834
## 10 59 125 1 0 0.89855072 0.02860313 0.9503632 0.84956295
## 11 60 124 1 0 0.89130435 0.02972717 0.9447781 0.84085714
## 12 65 123 2 0 0.87681159 0.03190746 0.9333961 0.82365740
## 13 71 121 1 0 0.86956522 0.03296902 0.9276100 0.81515252
## 14 81 120 1 0 0.86231884 0.03401448 0.9217668 0.80670491
## 15 88 119 2 0 0.84782609 0.03606427 0.9099232 0.78996675
## 16 92 117 1 0 0.84057971 0.03707173 0.9039292 0.78166991
## 17 93 116 1 0 0.83333333 0.03806935 0.8978906 0.77341763
## 18 95 115 1 0 0.82608696 0.03905833 0.8918099 0.76520757
## 19 105 114 1 0 0.81884058 0.04003974 0.8856890 0.75703763
## 20 107 113 1 0 0.81159420 0.04101457 0.8795299 0.74890594
## 省略一部分
## 192 550 15 1 0 0.34325958 0.18475887 0.4930486 0.23897674
## 193 551 14 0 1 0.34325958 0.18475887 0.4930486 0.23897674
## 194 559 13 0 1 0.34325958 0.18475887 0.4930486 0.23897674
## 195 588 12 0 1 0.34325958 0.18475887 0.4930486 0.23897674
## 196 641 11 1 0 0.31205416 0.20791044 0.4690333 0.20761384
## 197 654 10 1 0 0.28084875 0.23310483 0.4434980 0.17784977
## 198 687 9 1 0 0.24964333 0.26120251 0.4165393 0.14961805
## 199 705 8 1 0 0.21843791 0.29340057 0.3882139 0.12290937
## 200 728 7 1 0 0.18723250 0.33150176 0.3585552 0.09777018
## 201 731 6 1 0 0.15602708 0.37845310 0.3275969 0.07431220
## 202 735 5 1 0 0.12482167 0.43957565 0.2954319 0.05273786
## 203 740 4 0 1 0.12482167 0.43957565 0.2954319 0.05273786
## 204 765 3 1 0 0.08321444 0.59991117 0.2696771 0.02567754
## 205 821 2 0 1 0.08321444 0.59991117 0.2696771 0.02567754
## 206 965 1 0 1 0.08321444 0.59991117 0.2696771 0.02567754
## strata sex
## 1 sex=1 1
## 2 sex=1 1
## 3 sex=1 1
## 4 sex=1 1
## 5 sex=1 1
## 6 sex=1 1
## 7 sex=1 1
## 8 sex=1 1
## 9 sex=1 1
## 10 sex=1 1
## 11 sex=1 1
## 12 sex=1 1
## 13 sex=1 1
## 14 sex=1 1
## 15 sex=1 1
## 16 sex=1 1
## 17 sex=1 1
## 18 sex=1 1
## 19 sex=1 1
## 20 sex=1 1
## 省略一部分
## 192 sex=2 2
## 193 sex=2 2
## 194 sex=2 2
## 195 sex=2 2
## 196 sex=2 2
## 197 sex=2 2
## 198 sex=2 2
## 199 sex=2 2
## 200 sex=2 2
## 201 sex=2 2
## 202 sex=2 2
## 203 sex=2 2
## 204 sex=2 2
## 205 sex=2 2
## 206 sex=2 2
通过ggsurvplot()
进行可视化,非常多的细节可以修改,超级详细的教程可以参考后面的推文:超级详细的R语言生存分析可视化
ggsurvplot(fit.logrank,
data = df,
surv.median.line = "hv", # Add medians survival
# Change legends: title & labels
legend.title = "Sex",
legend.labs = c("Male", "Female"),
# Add p-value and tervals
pval = TRUE, # 这里P值直接写数字也行
conf.int = TRUE,
# Add risk table
risk.table = TRUE,
tables.height = 0.2,
tables.theme = theme_cleantable(),
ncensor.plot = TRUE,
# Color palettes. Use custom color: c("#E7B800", "#2E9FDF"),
# or brewer color (e.g.: "Dark2"), or ggsci color (e.g.: "jco")
palette = c("#E7B800", "#2E9FDF"),
ggtheme = theme_bw(), # Change ggplot2 theme
# Change font size, style and color
main = "Survival curve",
font.main = c(16, "bold", "darkblue"),
font.x = c(14, "bold.italic", "red"),
font.y = c(14, "bold.italic", "darkred"),
font.tickslab = c(12, "plain", "darkgreen")
)
自带的surv_cutpoint()
可用于寻找最佳切点,但是只能用于连续性数据。
使用myeloma
数据进行演示。
rm(list = ls())
# 0. Load some data
data(myeloma)
head(myeloma)
## molecular_group chr1q21_status treatment event time CCND1 CRIM1
## GSM50986 Cyclin D-1 3 copies TT2 0 69.24 9908.4 420.9
## GSM50988 Cyclin D-2 2 copies TT2 0 66.43 16698.8 52.0
## GSM50989 MMSET 2 copies TT2 0 66.50 294.5 617.9
## GSM50990 MMSET 3 copies TT2 1 42.67 241.9 11.9
## GSM50991 MAF <NA> TT2 0 65.00 472.6 38.8
## GSM50992 Hyperdiploid 2 copies TT2 0 65.20 664.1 16.9
## DEPDC1 IRF4 TP53 WHSC1
## GSM50986 523.5 16156.5 10.0 261.9
## GSM50988 21.1 16946.2 1056.9 363.8
## GSM50989 192.9 8903.9 1762.8 10042.9
## GSM50990 184.7 11894.7 946.8 4931.0
## GSM50991 212.0 7563.1 361.4 165.0
## GSM50992 341.6 16023.4 2096.3 569.2
寻找最佳切点:
# 1. Determine the optimal cutpoint of variables
res.cut <- surv_cutpoint(myeloma, time = "time", event = "event",
variables = c("DEPDC1", "WHSC1", "CRIM1") # 找这3个变量的最佳切点
)
summary(res.cut)
## cutpoint statistic
## DEPDC1 279.8 4.275452
## WHSC1 3205.6 3.361330
## CRIM1 82.3 1.968317
查看根据最佳切点进行分组后的数据分布情况:
# 2. Plot cutpoint for DEPDC1
plot(res.cut, "DEPDC1", palette = "npg")
## $DEPDC1
根据最佳切点重新划分数据,这样数据就根据最佳切点变成了高表达/低表达组。
# 3. Categorize variables
res.cat <- surv_categorize(res.cut)
head(res.cat)
## time event DEPDC1 WHSC1 CRIM1
## GSM50986 69.24 0 high low high
## GSM50988 66.43 0 low low low
## GSM50989 66.50 0 low high high
## GSM50990 42.67 1 low high low
## GSM50991 65.00 0 low low low
## GSM50992 65.20 0 high low low
根据最佳切点绘制生存曲线:
# 4. Fit survival curves and visualize
library("survival")
fit <- survfit(Surv(time, event) ~DEPDC1, data = res.cat)
ggsurvplot(fit, data = res.cat, risk.table = TRUE, conf.int = TRUE)
确定最佳切点的R包还有非常多,其他的后续会再介绍。
多时间点和多指标的ROC曲线绘制,可参考:R语言画多时间点ROC和多指标ROC曲线
平滑版ORC曲线和最佳截点:生存资料ROC曲线的最佳截点和平滑曲线
ROC曲线的比较:ROC(AUC)曲线的显著性检验
下次继续介绍Cox回归。
参考资料
http://www.sthda.com/english/wiki/survival-analysis-basics https://www.emilyzabor.com/tutorials/survival_analysis_in_r_tutorial.html survival包帮助文档 https://mp.weixin.qq.com/s/2rwxeaF_M0UnqPi2F9JNxA
🔖精选合集